﻿
using System;
using System.Collections.Generic;
using System.Linq;

namespace EmperialApps.WeatherSpark.Data {

    using Point = Pair<double, double>;

    /// <summary>A direct implementation of Akima interpolation.</summary>
    /// <see href="http://student.ndhu.edu.tw/~u9111023/akima.pdf"/>
    /// <see href="http://www.iue.tuwien.ac.at/phd/rottinger/node60.html"/>
    public sealed class AkimaInterpolation {

        /// <summary>Retrieves a scale for the specified data values, accounting for interpolation eccentricities.</summary>
        public static Scale GetInterpolatedScale( bool isInverted, double screenSize, DataValues values, double dataMargin, double floor = double.NegativeInfinity ) {
            return GetInterpolatedScale( isInverted, screenSize, values, values.Interval, dataMargin, floor );
        }

        /// <summary>Retrieves a scale for the specified data values, accounting for interpolation eccentricities.</summary>
        public static Scale GetInterpolatedScale( bool isInverted, double screenSize, IEnumerable<double> values, int interval, double dataMargin, double floor = double.NegativeInfinity ) {
            IEnumerable<double> scaleValues;
            if( interval == 1 )
                scaleValues = values;
            else {
                var interpolation = new AkimaInterpolation( values, interval );
                scaleValues = interpolation.Interpolate( ConvertValue.Identity ).Select( p => p.Item2 );
            }

            return new Scale( isInverted, screenSize, scaleValues, dataMargin, floor );
        }

        /// <summary>Initializes a new instance of the <see cref="AkimaInterpolation"/> class.</summary>
        public AkimaInterpolation( IEnumerable<double> data, int interval ) {
            this._interval = interval;

            // 1 ≤ i ≤ k
            var y = new OffsetArray( "y", 1, data );
            int k = y.Count;
            if( k == 0 ) {
                this._segments = new Segment[0];
                return;
            }

            // slopes
            var d = new OffsetArray( "d", -1, k + 4 );
            for( int i = 1; i < k; ++i )
                d[i] = y[i + 1] - Safe( y[i] );
            d[0] = 2 * d[1] - Safe( d[2] );
            d[-1] = 2 * d[0] - Safe( d[1] );
            d[k] = 2 * d[k - 1] - d[k - 2];
            d[k + 1] = 2 * d[k] - d[k - 1];
            d[k + 2] = 2 * d[k + 1] - d[k];

            // weights
            var w = new OffsetArray( "w", -1, k + 4 );
            for( int i = 0; i <= k + 2; ++i )
                w[i] = Math.Abs( Safe( d[i] ) - Safe( d[i - 1] ) );

            // derivative estimates
            var t = new OffsetArray( "t", 0, k + 2 );
            for( int i = 1; i <= k + 1; ++i ) {
                double d1 = d[i - 1];
                double d2 = d[i];
                bool previousClose = AreClose( d1, d[i - 2] );
                bool nextClose = AreClose( d2, d[i + 1] );

                double v;
                if( nextClose && !previousClose )
                    v = d2;
                else if( (previousClose && !nextClose) || AreClose( d1, d2 ) )
                    v = d1;
                else if( previousClose && nextClose )
                    v = (d1 + d2) / 2;
                else {
                    double w1 = w[i + 1];
                    double w2 = w[i - 1];
                    v = (d1 * w1 + d2 * w2) / (w1 + w2);
                }

                t[i] = Safe( v );
            }

            // interpolation segments
            this._segments = new Segment[k];
            for( int i = 1; i <= k; ++i )
                this._segments[i - 1] = new Segment( y[i], d[i], t[i], t[i + 1] );
        }

        /// <summary>Gets the number of source data values used to create the interpolation.</summary>
        public int Count { get { return this._segments.Length; } }

        /// <summary>Gets the source data value at the specified index.</summary>
        public double this[int index] { get { return this._segments[index].Sample; } }

        /// <summary>Returns the interpolated data and range values between the specified range points.</summary>
        public IEnumerable<Point> Interpolate( int index, double x1, double x2 ) {
            Segment segment = this._segments[index];
            double start = Math.Ceiling( x1 );
            double end = Math.Ceiling( x2 );
            double range = end - start;

#if DEBUG
            double offset = 1;
            double y1 = segment.Interpolate( 0 );
            if( !double.IsNaN( y1 ) && this._interval > 1 ) {
                yield return new Point( x1 - offset, y1 + offset );
                yield return new Point( x1 + offset, y1 + offset );
                yield return new Point( x1 + offset, y1 - offset );
                yield return new Point( x1 - offset, y1 - offset );
                yield return new Point( x1 - offset, y1 );
            }
#endif
            for( int i = 0; i < range; i += Spacing ) {
                double x = start + i;
                double y = segment.Interpolate( i / range );
                if( !double.IsNaN( y ) )
                    yield return new Point( x, y );
            }
        }

        /// <summary>Returns the interpolated data and range values for all of the data points used to create the interpolation.</summary>
        public IEnumerable<Point> Interpolate( Scale scale, int offset = 0 ) {
            return Interpolate( scale.ToScreen, offset );
        }

        /// <summary>Returns the interpolated data and range values for all of the data points used to create the interpolation.</summary>
        public IEnumerable<Point> Interpolate( Func<double, double> toScreen, int offset = 0 ) {
            int count = Math.Max( 0, this.Count - 1 );
            return
                from i in Enumerable.Range( 0, count )
                let x1 = toScreen( i * this._interval + offset )
                let x2 = toScreen( (i + 1) * this._interval + offset )
                from point in this.Interpolate( i, x1, x2 )
                select point;
        }

        #region Private Members

        private const int Spacing = 4;

        private readonly Segment[] _segments;
        private readonly int _interval;

        private static double Safe( double value ) {
            return double.IsNaN( value ) ? 0 : value;
        }

        private static bool AreClose( double a, double b ) {
            return Math.Abs( a - b ) < 1e-6;
        }

        private struct Segment {
            private readonly double[] _coefficients;

            public Segment( double sample, double d1, double t1, double t2 ) {
                double c1 = t1;
                double c2 = 3 * d1 - 2 * t1 - t2;
                double c3 = t1 + t2 - 2 * d1;

                if( AreClose( c3, 0 ) && AreClose( c2, 0 ) )
                    this._coefficients = new[] { sample, c1 };
                else if( AreClose( c3, 0 ) )
                    this._coefficients = new[] { sample, c1, c2 };
                else
                    this._coefficients = new[] { sample, c1, c2, c3 };
            }

            public double Sample { get { return this._coefficients[0]; } }

            public double Interpolate( double delta ) {
                // s(x) = c₀ + c₁(x - xᵢ) + c₂(x - xᵢ)² + c₃(x - xᵢ)³
                // delta = (x - xᵢ)
                double result = 0;
                double product = 1;
                foreach( double c in this._coefficients ) {
                    result += c * product;
                    product *= delta;
                }

                return result;
            }
#if DEBUG
            public override string ToString( ) {
                return "Segment " + this.Sample;
            }
#endif
        }

        private struct OffsetArray {
            private readonly string _name;
            private readonly int _offset;
            private readonly double[] _array;

            public OffsetArray( string name, int offset, IEnumerable<double> values ) {
                this._name = name;
                this._offset = offset;
                this._array = values.ToArray( );
            }

            public OffsetArray( string name, int offset, int length ) {
                this._name = name;
                this._offset = offset;
                this._array = new double[length];
            }

            public int Count { get { return this._array.Length; } }

            public double this[int index] {
                get { return this._array[this.GetArrayIndex( index )]; }
                set { this._array[this.GetArrayIndex( index )] = value; }
            }

            public override string ToString( ) {
                return this._name + ": Count=" + this._array.Length + ", Offset=" + this._offset;
            }

            private int GetArrayIndex( int index ) {
                int arrayIndex = index - this._offset;
                if( (uint)arrayIndex >= this._array.Length )
                    throw new IndexOutOfRangeException( string.Format( "Index {0} was outside the range [{1},{2}] of '{3}'.", index, this._offset, this.Count - 1 + this._offset, this._name ) );

                return arrayIndex;
            }
        }

        #endregion

    }

}
